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HIGHLIGHTS 


•  A  PEMFC  cold  start  stack  model  is  developed. 

•  A  novel  variable  heating  and  load  control  (VHLC)  method  is  proposed. 

•  Dead  individual  cells  can  be  activated  again  by  other  cells  with  VHLC. 

•  Proper  VHLC  improves  cold  start  performance  significantly. 
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In  this  study,  a  cold  start  model  for  proton  exchange  membrane  fuel  cell  (PEMFC)  stacks  is  developed, 
and  a  novel  start-up  method,  variable  heating  and  load  control  (VHLC),  is  proposed  and  evaluated.  The 
main  idea  is  to  only  apply  load  to  the  neighboring  still-active  cells,  and  to  apply  external  heating  to 
certain  cells  inside  the  stack  simultaneously  (load  is  not  applied  to  the  cells  fully  blocked  by  ice,  although 
these  cells  can  gain  heat  from  neighboring  cells).  With  the  VHLC  method,  it  is  found  that  the  stack 
voltage  first  increases,  then  decreases  due  to  the  full  blockage  of  ice  in  some  of  the  individual  cells,  and 
finally  the  dead  cells  are  heated  by  the  other  active  cells  and  activated  again  one  by  one.  Based  on  this 
method,  the  external  heating  power  and  the  stack  self-heating  ability  are  utilized  more  efficiently.  With 
proper  implementation  of  the  VHLC  method,  it  is  demonstrated  that  the  cold  stat  performance  can  be 
improved  significantly,  which  is  critically  important  for  PEMFC  in  automotive  applications. 

©  2013  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Although  proton  exchange  membrane  fuel  cell  (PEMFC)  is 
considered  as  one  of  the  most  suitable  energy  conservation  devices 
in  the  future  for  automotive  applications,  there  are  still  some 
technical  issues  to  be  solved  before  its  commercial  application. 
Among  the  different  technical  issues,  start-up  from  subzero  tem¬ 
peratures,  known  as  “cold  start”,  is  an  important  one,  because  the 
ice  formed  during  the  cold  start  processes  can  block  the  reaction 
sites,  hindering  the  reactant  transport  and  causing  material 
degradation. 

Most  of  the  previous  experimental  studies  of  PEMFC  cold  start 
were  conducted  for  performance  degradation  measurement  [1—6], 
ice  formation  visualization  [7-9]  and  property  characterization 
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[10-15].  On  the  other  hand,  mathematical  models  are  needed  to 
study  the  detailed  transport  processes.  One-dimensional  models 
and  analytical  models  were  first  developed  to  study  the  variations 
of  voltage,  current  and  temperature  during  the  cold  start  processes 
[16-18],  Sundaresan  and  Moore  [19]  developed  a  one-dimensional 
thermal  model  for  simulating  the  cold  start  processes  of  PEMFC 
stack,  and  this  model  separated  the  fuel  cell  into  different  layers  to 
obtain  the  one-dimensional  temperature  distribution  in  a  stack. 
Another  one-dimensional  cold  start  model  for  PEMFC  stack  was 
developed  by  Khandelwal  et  al.  [20]  to  investigate  the  effects  of 
various  operating  conditions  and  external  heating  methods  on  the 
start-up  performance.  Moreover,  more  detailed  transport  phe¬ 
nomena  during  the  cold  start  processes  of  PEMFC  were  investigated 
based  on  multi-dimensional  multiphase  models  in  recent 
years  [21-28],  For  example,  Jiao  and  Li  [25,26]  developed  a 
three-dimensional  multiphase  model  to  simulate  the  cold  start 
processes  in  a  single  PEMFC  with  straight  flow  channels,  and 
investigated  the  various  transport  characteristics.  Guo  et  al.  [27] 
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Nomenclature 

BP 

bipolar  plate 

c 

cathode 

A 

cell  geometric  area,  m2 

cell 

cell  characteristic 

ASR 

area  specific  resistance,  fi  m2 

channel 

flow  channel 

c 

mole  concentration,  mol  m-3 

CL 

catalyst  layer 

Cp 

specific  heat,  J  kg  1  l<  1 

cone 

concentration 

D 

mass  diffusivity,  m2  s_1 

Cond 

condensation 

EW 

equivalent  weight  of  membrane,  kg  kmol  1 

eff 

effective 

F 

Faraday’s  constant,  96487  C  mol1 

eq 

equilibrium 

h 

surrounding  heat  transfer  coefficient,  W  m  2  1<  1 

/ 

frozen 

I 

current  density,  A  cm-2 

fl 

fluid  phase 

j 

reaction  rate,  A  m  3 

fmw 

frozen  membrane  water 

k 

thermal  conductivity,  W  m  1  K_1 

GDL 

gas  diffusion  layer 

M 

molecular  weight,  kg  mol1 

h2 

hydrogen 

P 

pressure,  Pa 

h2o 

water 

Q 

heat  transfer  rate,  W 

i.j 

the  ith  and  jth  components 

R 

universal  gas  constant,  J  mol-1  K_1 

ice 

ice 

S 

volume  fraction 

Kns 

Knudsen  diffusion  coefficient 

S 

source  terms 

iq 

liquid  water 

t 

time,  s 

l-i 

liquid  water  to  ice  (vice  versa) 

T 

temperature,  K 

mem 

membrane 

To 

volume  averaged  cell  temperature,  K 

nemst 

Nemst 

nf 

non-frozen 

Creek  letters 

nmw 

non-frozen  membrane  water 

a 

transfer  coefficient 

n-f 

non-frozen  membrane  water  to  frozen  membrane 

e 

porosity 

water  (vice  versa) 

K 

water  transfer  rate,  s_1 

n—i 

non-frozen  membrane  water  to  ice  (vice  versa) 

V 

overpotential,  V 

n-v 

non-frozen  membrane  water  to  vapor  (vice  versa) 

A 

water  content  in  ionomer 

o2 

oxygen 

? 

stoichiometry  ratio 

ohmic 

ohmic 

P 

density,  kg  m  1 

out 

outlet 

(0 

volume  fraction  of  ionomer  in  catalyst  layer 

pc 

phase  change 

8 

thickness,  m 

r 

mean  pore  radius 

ref 

reference  state 

Subscripts  and  superscripts 

sat 

saturation 

a 

anode 

si 

solid  phase 

act 

activation 

T 

energy  (for  source  term) 

amb 

ambient 

vp 

water  vapor 

atm 

atmosphere 

wall 

surrounding  wall  of  the  stack 

further  implemented  the  assisted  catalytic  hydrogen-oxygen  reac¬ 
tion  method  in  the  multiphase  cold  start  model,  and  found  that 
successful  start-up  from  -20  °C  can  be  achieved  with  proper 
control  of  this  method. 

It  should  be  noticed  that  even  though  the  multi-dimensional 
multiphase  cold  start  models  for  PEMFC  can  be  used  to  investi¬ 
gate  the  various  detailed  transport  phenomena,  these  models  were 
only  developed  for  single  PEMFCs  with  single  straight  flow  chan¬ 
nels  (only  a  small  part  of  a  practical  single  PEMFC)  [21-27]  or  small 
stacks  up  to  only  3  single  cells  [28]  to  the  best  of  the  authors’ 
knowledge,  although  previous  studies  showed  that  the  transport 
characteristics  might  be  different  at  larger  scales  [19,20,28—31], 
The  main  restriction  of  multi-dimensional  multiphase  cold  start 
models  for  PEMFC  is  the  limitation  of  computational  power  [28],  On 
the  other  hand,  as  mentioned  previously,  although  one¬ 
dimensional  models,  requiring  insignificant  computational  power, 
might  be  an  alternative  way  to  study  the  cold  start  characteristics  of 
PEMFC  stacks,  these  models  often  focused  on  the  heat  generation 
and  transfer  processes  [19,20],  and  the  coupled  heat  and  mass 
transfer  (e.g.  ice  formation  and  water  transport)  was  often  simpli¬ 
fied.  The  major  limitation  of  one-dimensional  model  is  that  the 
detailed  transport  phenomenon  in  the  other  two  directions  cannot 


be  taken  into  account,  for  example,  the  effect  of  detailed  flow  field 
structure  [32]  and  water  distributions  cannot  be  considered.  Higher 
dimensional  models  are  needed  to  satisfy  such  requirements,  and 
the  results  of  previous  multi-dimensional  cold  start  models  figured 
out  these  detailed  transport  phenomenons  for  single  PEMFCs  [21— 
28],  However,  for  PEMFC  stacks,  considering  the  computational 
power  limitation,  and  the  fact  that  the  through-plane  transport 
plays  the  leading  role,  one-dimensional  model  is  a  reasonable 
choice  for  PEMFC  stack  modeling.  Therefore,  it  is  necessary  to 
develop  a  one-dimensional  (to  be  computationally  efficient)  large- 
scale  cold  start  model  considering  the  major  heat  and  mass  transfer 
processes  for  PEMFC  stacks.  With  such  a  comprehensive  stack 
model,  the  various  assisted  start-up  methods  can  also  be  evaluated 
and  analyzed  in  details. 

In  this  study,  a  one-dimensional  cold  start  model  is  developed 
for  PEMFC  stacks.  This  model  accounts  for  the  major  transport 
phenomena  during  both  the  failed  and  successful  cold  start  pro¬ 
cesses,  such  as  the  coupled  water  transport/phase  change  and  heat 
transfer.  Based  on  this  model,  a  novel  start-up  method,  variable 
heating  and  load  control  (VHLC),  is  proposed  and  evaluated.  The 
transport  characteristics  during  the  different  cold  start  processes  in 
different  PEMFC  stacks  are  investigated. 
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2.  Mathematical  model 

Fig.  1  shows  the  major  transport  phenomena  considered  in  this 
stack  model,  including  the  water  transport  and  phase  change  in 
membrane,  catalyst  layer  (CL)  and  gas  diffusion  layer  (GDL),  and 
heat  generation/transfer  in  different  components.  The  heat  and 
mass  transfer  along  the  direction  normal  to  the  membrane  surface 
(the  through-plane  direction)  is  considered  in  this  one¬ 
dimensional  model  (defined  as  the  x-direction  in  this  model).  The 
design  and  operating  parameters  are  given  in  Table  1. 

2.1.  Water  transport 

In  this  model,  for  the  water  transport  in  the  membrane  and  in 
the  electrolyte  of  CL,  the  mass  conservation  equations  of  non- 
frozen  and  frozen  membrane  water  are  solved  along  the  through- 
plane  direction  (x-direction): 


(Table  2).  Snmw  and  Sfmw  are  the  source  terms  of  non-frozen 
membrane  water  and  frozen  membrane  water,  respectively,  as 
shown  in  Table  3.  The  source  term  of  non-frozen  membrane 
water  includes  the  product  water  from  the  reaction,  water  phase 
change,  and  migrated  water  due  to  the  electro-osmotic  drag 
(EOD)  effect.  The  source  term  of  frozen  membrane  water  in¬ 
cludes  the  phase  change  between  frozen  and  non-frozen  mem¬ 
brane  water.  These  assumptions  of  water  phase  change  were 
proposed  in  previous  work  for  single  PEMFC  modeling  [25],  In 
the  pores  of  CL  and  GDL,  the  mass  conservation  equations  of 
liquid  water  and  ice  are  solved  along  the  through-plane  direction 
(x-direction): 

0(ePlq'slq)  92  (PlqD|q'slq) 

Liquid  water  :  ^ — L  =  v  ^ - L  +  slq  (3) 


Non  -  frozen  membrane  water  : 


®(w^nf) 


Pmem  ^  (^^nmw^nf) 

EW  9X2 


(1) 


Frozen  membrane  water  :  ■ 


(4) 


where  p  (kg  m-3)  is  the  density,  EW  (kg  kmol-1)  the  equivalent 
weight  of  membrane,  w  the  volume  fraction  of  ionomer  in  CL,  X 
the  water  content  in  ionomer,  and  D  (m2  s_1)  the  mass  diffusivity 


where  e  is  the  porosity,  and  s  the  volume  fraction.  Since  the 
convective  mass  transfer  in  GDL  and  CL  is  only  considerable  with 
cross-flow  along  the  in-plane  directions  [33],  only  diffusive  mass 


Fig.  1.  Major  transport  phenomena  considered  in  this 


model. 
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Kcl  =  6.2  x  10“13  m2;  KGDl  =  6.2  x  10  12  m 


5m  =  2.0 

C  =  W  =  -20  °C 
Unit  =  rsurr  =  -20  °C 
P™'  =  1  atm 


h  =  2  W  m“2  K  1 
0 

2.5 


h  =  0  W  nr2  K-1 


transfer  is  considered  in  this  one-dimensional  model  along  the 
through-plane  direction.  The  related  transport  properties  are  given 
in  Table  2.  The  source  terms  consider  the  water  phase  change,  as 
shown  in  Table  3.  The  water  vapor  transport  is  also  considered  to  be 
the  diffusion  dominated  along  the  through-plane  direction,  and  the 
average  water  vapor  concentration  in  CL  and  GDL  is  updated 
directly  at  each  time  step  rather  than  solving  a  partial  differential 
equation.  In  fact,  during  the  cold  start  processes,  the  water  vapor 
amount  is  very  low  due  to  the  low  water  saturation  pressure.  Such 
low  water  vapor  amount  has  insignificant  effect  on  the  cold  start 
processes  (this  was  also  reported  in  previous  studies  [25,26,33]). 
Therefore,  the  water  vapor  transport  is  further  simplified  to 
improve  the  calculation  efficiency.  For  water  vapor  in  CL: 


Cvp,CL  —  Cvp^CL  +  ? (2'-cl  ~~  i© 

(4s"-4cl,r)D«CD1  At 
la.*®. 


where  c  (mol  m-3)  is  the  concentration,  t  (s)  the  time  step,  At  (s)  the 
time  step  size,  C  (s_1)  the  water  transfer  rate,  <5  (m)  the  thickness,  ?  the 
stoichiometry  ratio,  and  D  (m2  s_1)  the  effective  diffusion  coefficient 
corrected  by  the  ice  and  liquid  volume  fraction  in  CL  and  GDL  (Table  2). 
The  second  term  on  the  right  hand  side  of  Equation  ( 5 )  represents  the 
water  absorption/desorption  by  the  ionomer  in  CL  and  the  third  term 
is  for  water  diffusion  between  CL  and  GDL  For  water  vapor  in  GDL: 


Parameters 

Value 

Membrane  water  diffusivity  (m2  s-1) 

'  2.69266  xlO-10(7<  2) 

10-1°  exp [241 6 (3 J3  - 1)]  ■  [0.87.  (3  -  7nf)  +  2.95- (fnf  -  2)] ,  (2  <  7^  <  3) 

10-1°. exp[2416(3J3 -].)].  [2.95-(4-Anf)  +  1.642454- (7nt -3}J,  (3  <  ^  <  4) 

10-10 •  exp [241 6 (353  -  g] ■  [2.563  -  0.33 -7nf  +  0.0264-4  -  0.000671  -4)],  (4  >  4) 

Liquid  water  diffusivity  (m2  s-1) 

Dlq  = 

^=Kc°L,GD^q(l-%e)4 

Water  vapor  diffusivity  in  CL  (m2  s-1) 

D,  =  2.982  xlO-^)^(Mp) 

Water  vapor  diffusivity  in  GDL  (mz  s-1) 

DfDL=D^(l-sPq--sP-)''5 
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er  conservation  equations. 


Non-frozen  membrane 
(kmol  m  5  s1) 

Frozen  membrane 


-Sn_f  (in  membrane) 

-  Sn-v  -  S„_ ,*  +  SE0D  (in  cathode  CL) 
-S„-„  -  S„_i  +  Shod  (in  anode  CL) 


rate  is  used  in  this  model,  the  reaction  rate  is  calculated  as  the  ratio 
of  the  current  density  (/,  A  m-2)  and  CL  thickness  (<5cl.  na).  The 
relationship  between  the  current  density  and  stack  voltage  is  dis¬ 
cussed  in  the  next  subsection. 


2.3.  Performance  prediction 

In  this  model,  the  charge  conservation  equation  is  not  solved, 
and  as  the  alternative,  the  performance  is  calculated  analytically 
based  on  the  Tafel  representation  [34]  and  the  transport  parame¬ 
ters  during  the  cold  start  processes.  The  output  voltage  can  be 
calculated  as 


Vout  —  Vnernst  +  Vact  +  Vconc  +  ^ohmic 


(9) 


A  _  t_At  (CCLAt  CGDLf)  ^CL-GDL  At 

v  P'GDL  =  Cvp’GDL  +  ^S72)  +  (G72n^LCCL 

(CGDLf  —  cchannel)DGDL  M  ^ 

^gdl/2  <5gdl«gdl 

where  the  second  and  third  terms  on  the  right  hand  side  of  Equation 
(6)  represent  the  water  diffusion  between  CL  and  GDL  and  between 
GDL  and  flow  channel,  respectively.  As  mentioned  previously,  since 
the  water  removed  as  vapor  counts  only  a  very  small  proportion  of 
product  water  (this  was  also  reported  in  previous  studies  [25,26,33]), 
it  is  assumed  that  the  water  vapor  is  taken  away  instantaneously  in 
the  flow  channel,  therefore  the  water  vapor  concentration  in  flow 
channel  (cchannei)  is  considered  to  be  zero  in  this  model. 


2.2.  Heat  transfer 

The  one-dimensional  heat  transfer  along  the  through-plane 
direction  (x-direction)  with  different  heating  sources  is  consid¬ 
ered  in  this  model: 


d2( keff  t) 

Energy  conservation:  ^((pCp^r)  =  fo2  +Sr  (7) 

where  Cp  (J  kg-1  K-1)  is  the  specific  heat,  T  (K)  the  temperature,  and 
k  (W  m  1  K-1)  the  thermal  conductivity.  The  convective  heat 
transfer  along  the  through-plane  direction  is  neglected,  since  it  is 
often  insignificant  by  comparing  with  the  conductive  heat  transfer. 
The  effective  density,  heat  capacity  and  thermal  conductivity  are 
calculated  based  on  the  stack  materials  and  ice/liquid/gas  concen¬ 
trations.  The  heating  sources  in  the  different  components  include 
the  activational  heat,  ohmic  heat,  reversible  heat  and  latent  heat,  as 
demonstrated  in  Fig.  1.  The  source  term  is  expressed  as 


'jcV act+^^  +  V  (in  anode  CL) 

— Jc2f S  +jcV act  + 1  3^C1‘  +  Spc  (in  cathode  CL) 
^r+5pc  (in  GDL) 

(in  BP) 

;2ASRmem  +  spc  (in  membrane) 

,  0  (in  other  zones) 


(8) 


where  ASR  (Q  m2)  is  the  area  specific  resistance,  t]  (V)  the  over¬ 
potential,  and  j  (A  m-3)  the  reaction  rate.  Since  the  average  reaction 


The  first  term  on  the  right  hand  side  of  Equation  (9)  is  the  open 
circuit  voltage  (OCV),  which  is  calculated  by  the  Nernst  equation: 

Vnernst  » 1-23  -  0.9  x  KT3(T  -  T0)  +  (10) 

2c  V  Ph20  / 


where  R  (J  mol1  K_1)  is  the  universal  gas  constant,  and  T0  is  a 
constant  of  298.15  K.  The  second  term  on  the  right  hand  side  of 
Equation  (9)  is  the  activational  voltage  loss: 


Vact  =  -b  In ) 


<l-sic 


b 


RT 

aF 


(11) 


(12) 


(13) 


1+  (1-1/?)  0.21  pc 

2  RT 


(14) 


where  b  is  the  Tafel  slope,  /  (A  cm-2)  is  the  current  density,  j5 
(A  cm-2)  is  the  reference  current  density.  Therefore,  the  effects  of 
reactant  concentration  and  ice/liquid  formation  on  the  activational 
voltage  loss  are  all  considered.  The  third  term  on  the  right  hand  side 
of  Equation  (9)  is  the  voltage  loss  due  to  the  reactant  transport: 

(15) 


Jd 


4 Fch 

SGDL/DfDl  +  0.55cl/D^ 


(16) 


The  fourth  term  on  the  right  hand  side  of  Equation  (9)  is  the 
ohmic  voltage  loss: 

^ohmic  =  -ASR  j  (17) 

where  the  area  specific  resistance  (ASR)  is  calculated  as 

ASR  =  Ace||(Rmem  +  Rcl  +  Pgdl  +  rbp)  (18) 


The  resistance  of  membrane  is  calculated  based  on  the  corre¬ 
lation  provided  by  Springer  et  al.  [35],  and  the  resistance  of  CL  in¬ 
cludes  both  the  ionomer  (for  proton  transport)  and  catalyst  (for 
electron  transport). 
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2.4.  Boundary  and  initial  conditions 

For  the  surrounding  walls  of  the  stack,  the  convective  heat 
transfer  condition  is  set  as: 

(L  =  M(ramb  -  rwaU)  (19) 

where  h  (W  m  2  I<)  is  the  heat  transfer  coefficient,  A  (m2)  the 
surface  area,  Tamb  (K)  the  ambient  temperature,  and  Twau  (K)  the 
wall  temperature.  The  heat  transfer  coefficient  is  set  to  be 
2  W  m-2  K  on  the  side  walls  except  the  end  plates,  to  reflect  a  well 
insulation  condition  on  these  walls  (nature  convection  condition). 
On  the  end  plates,  different  heat  transfer  coefficients  are  consid¬ 
ered,  as  given  in  Table  1.  External  heating  is  also  applied  to  different 
cells  for  the  assisted  cold  start  processes. 

To  calculate  the  cell  output  voltage  at  a  current  density  (con¬ 
stant  current  condition),  the  current  density  of  each  cell  is 
considered  as  a  constant: 

fceii  j  =  Constant  (20) 

where  i  represents  the  ith  cell.  The  total  stack  output  voltage  is  the 
sum  of  all  cell  voltages  (can  be  calculated  by  using  Equation  (9)).  For 
the  constant  voltage  condition,  the  total  voltage  of  the  different 
cells  is  a  constant: 

E^uti  =  Wac  Inconstant  (21) 

For  a  stack  with  N  cells,  the  following  equation  system  is  solved: 


Vout,  -* 

FunTafel 

(I,  [ Parameters i]) 

Vout,  = 

FunTafel 

(/,  [ParameterSi]) 

^outN  =* 

PunTafel 

(/,  [[ ParametersN ]) 

EVi  = 

Vout_stack 

where  Funjafei  is  the  performance  equation  (Equation  (9)),  and  / 
(A  m-2)  is  the  current  density  (one  of  the  unknowns).  ParameterSi  is 
the  state  parameters  of  the  ith  individual  cell  such  as  temperature 
and  ice  volume  fraction.  Vout ,  (V)  is  the  output  cell  voltage  of  the  ith 
cell  (one  of  the  unknowns).  Therefore,  the  current  density  and  the 
voltages  of  different  cells  can  be  calculated. 

The  initial  temperature  is  set  as  the  same  with  the  ambient 
temperature  of  -20  °C,  the  initial  water  content  in  ionomer  is  2.5  to 
reflect  the  well  purged  condition  before  the  cold  start  processes, 
and  the  initial  ice  volume  fraction  is  0. 


2.5.  Numerical  methods 

In  this  study,  the  membrane  water,  liquid,  ice  and  energy  con¬ 
servation  equations  are  solved.  The  finite  volume  method  is  used  to 
discretize  these  differential  equations.  There  are  totally  250  grid 
nodes  for  each  cell,  and  the  distributions  for  the  different  compo¬ 
nents  are  30, 20, 20, 30  and  50  nodes  for  the  bipolar  plate  (BP),  flow 
channel,  GDL,  CL  and  membrane,  respectively.  In  each  discretized 
control  volume,  the  equations  are  solved  with  the  fourth-order 
Runge-Kutta  method.  Explicit  time  discretization  is  used,  and 
adaptive  time  step  using  the  Dormand-Prince  pair  for  error  esti¬ 
mation  is  adopted  for  the  automatic  adjustment  of  iteration  time 
step.  Grid  independent  study  has  been  carefully  carried  out  to 
ensure  that  further  increasing  the  number  of  grid  has  negligible 


influence  on  the  results.  To  solve  Equation  (22)  for  the  constant 
voltage  condition,  the  algorithm  of  Levenberg— Marquardt  [36]  is 
used  with  the  second  order  accuracy.  The  termination  tolerance  is 
set  as  1(T9,  which  means  that  the  convergence  completes  once  the 
errors  become  less  than  10  9. 

3.  Results  and  discussion 

Fig.  2  shows  the  comparison  between  the  present  model  pre¬ 
diction  and  experimental  data  in  Ref.  [37]  for  the  cold  start  pro¬ 
cesses  starting  from  -20  °C  at  3  different  current  densities.  For  this 
comparison,  the  model  parameters  are  set  according  to  the 
experimental  conditions  in  Ref.  [37],  It  can  be  noticed  that  the 
agreement  is  reasonable  for  the  whole  cold  start  processes. 

In  this  section,  the  unassisted  cold  start  characteristics  of 
different  PEMFC  stacks  are  discussed  first,  and  the  effects  of  num¬ 
ber  of  cells  in  stacks  and  start-up  modes  (constant  current  and 
constant  voltage)  are  investigated.  The  external  heating  assisted 
cold  start  processes  are  then  analyzed,  and  the  external  heating 
method  is  optimized  to  improve  the  cold  start  performance.  A 
novel  start-up  method,  variable  load  and  heating  control  (VLHC),  is 
proposed  and  evaluated  as  well.  As  shown  in  Table  1,  a  typical  start¬ 
up  temperature,  -20  °C,  is  considered  for  all  the  cold  start  pro¬ 
cesses:  the  inlet  gases  are  also  at  -20  °C  and  the  relative  humidity 
is  0;  the  initial  water  content  in  electrolyte  is  2.5  and  there  is  no  ice 
at  the  beginning  of  the  cold  start  processes;  the  start-up  current 
density/voltage  are  0.1  A  cm-2,  0.15  A  cm^2  and  0.3  V;  and  the 
operating  pressure  is  1  atm.  At  the  two  end  surfaces,  forced  con¬ 
vection  with  a  heat  transfer  coefficient  of  100  W  m  2  K-1  is 
considered,  and  at  the  other  outer  surfaces,  a  heat  transfer  coeffi¬ 
cient  of  2  W  m-2  K-1  is  considered  representing  the  nature  con¬ 
vection  condition.  The  heat  transfer  coefficient  can  also  be  changed 
to  0  representing  the  heat  insulation  condition.  For  the  external 
heating  condition,  a  heating  power  of  100  W  is  applied  for  different 
cells  in  the  20-cell  stack  to  optimize  the  cold  start  performance. 

3.1.  Effect  of  number  of  cells 

Fig.  3  shows  the  effect  of  stack  cell  numbers  (total  amount  of 
cells  in  a  stack)  on  the  time  durations  and  final  volume  averaged 
stack  temperatures  for  the  failed  unassisted  cold  start  processes 
from  -20  °C  at  0.1  A  cm-2.  For  the  stacks  with  more  cells,  the  start¬ 
up  process  can  last  longer  and  the  final  temperature  reached  is 


Time,  s 


Fig.  2.  Comparison  between  the  present  model  predictions  and  experimental  data 
(experimental  data  from  Ref.  [37]). 
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Fig.  3.  Effect  of  number  of  cells  (total  amount  of  cells  in  a  stack)  on  the  time  durations 
and  final  volume  averaged  stack  temperatures  for  the  failed  unassisted  cold  start 
processes  from  -20  °C  at  0.1  A  cnr2. 


higher,  because  a  stack  with  more  cells  reduces  the  outer  surface 
area,  resulting  in  better  heat  insulation.  However,  this  effect  be¬ 
comes  less  significant  by  further  increasing  the  amount  of  cells.  The 
start-up  durations  and  final  temperatures  are  similar  for  the  stacks 
with  more  than  20  cells. 

The  temperature  distributions  in  the  stacks  with  different 
amounts  of  cells  (stack  with  1  cell  represents  single  cell)  at  the  end 
in  of  the  failed  unassisted  cold  start  processes  from  -20  °C  at 
0.1  A  cm  2  are  shown  in  Fig.  4.  For  the  different  stacks,  the  cells  in 
the  middle  have  the  highest  temperature,  on  the  other  hand,  the 
peak  temperatures  can  be  observed  around  the  membrane  in 
different  cells,  because  most  of  the  heat  is  produced  there.  The 
temperature  differences  between  the  cells  on  the  side  and  in  the 
middle  are  larger  for  the  stacks  with  more  cells. 

For  the  20-cell  stack  starting  from  -20  °C  at  0.1  A  cm  2,  the  ice 
volume  fractions  in  the  cathode  CLs  of  different  cells  at  the  end  of 
the  failed  unassisted  cold  start  process  is  shown  in  Fig.  5.  The  cells 
at  the  two  ends  have  slightly  more  ice  than  the  middle  cells,  mainly 
due  to  the  lower  temperatures  (Fig.  4)  there.  The  10  cells  in  the 
middle  have  almost  the  same  amount  of  ice.  Generally,  the  ice 
amounts  in  the  different  cells  of  the  stack  are  similar. 


Fig.  4.  Temperature  distributions  in  the  stacks  with  different  amounts  of  cells  (stack 
with  1  cell  represents  single  cell)  at  the  end  in  of  the  failed  unassisted  cold  start 
processes  from  -20  °C  at  0.1  A  cnrr2. 


Cell  Number 


Fig.  5.  Ice  volume  fractions  in  the  cathode  CLs  of  different  cells  in  the  20-cell  stack  at 
the  end  of  the  failed  unassisted  cold  start  process  from  -20  °C  at  0.1  A  cm-2. 


3.2.  Effect  of  start-up  current/voltage 

Fig.  6  shows  the  evolutions  of  average  cell  voltage/current 
density  and  average  ice  volume  fraction  in  cathode  CLs  of  different 
stacks  starting  from  -20  °C.  For  the  constant  current  start-up  mode 
at  0.15  A  cm-2  (Fig.  6(a)),  the  voltage  is  the  average  value  of  the 
different  cells  in  the  stack.  The  voltages  are  initially  almost  the 
same  for  the  different  stacks  with  different  amounts  of  cells.  After 
about  10  s,  the  difference  in  voltage  becomes  noticeable,  and  the 
stacks  with  more  cells  have  higher  voltages,  due  to  the  higher 
temperatures  caused  by  the  better  heat  insulation,  as  discussed 
previously  with  Figs.  3  and  4.  Moreover,  the  voltage  difference  is 
larger  between  the  5-cell  stack  and  single  cell  than  between  the 
100-cell  stack  and  5-cell  stack,  following  the  same  trend  of  the 
stack  temperature  shown  in  Fig.  3.  The  ice  volume  fractions  are 
quite  similar  for  the  different  stacks  during  the  whole  process,  and 
the  stacks  with  more  cells  only  have  slightly  slower  ice  formation 
processes.  For  the  constant  voltage  mode  (Fig.  6(b)),  it  can  be 
noticed  that  for  the  stacks  with  slightly  lower  ice  formation  rates 
(smaller  stacks),  the  current  densities  are  lower,  and  correspond¬ 
ingly  the  start-up  processes  last  slightly  longer.  These  trends  are 
opposite  to  Fig.  6(a)  with  constant  current  density.  Because  the 
larger  stacks  in  Fig.  6(b)  (higher  average  temperatures  as  discussed 
previously)  lead  to  higher  current  densities  during  the  start-up 
processes,  causing  more  water  generation  and  ice  formation. 

The  ice  volume  fractions  in  the  cathode  CLs  of  different  cells  in 
the  50-cell  stack  during  the  failed  unassisted  cold  start  processes 
from  -20  °C  at  0.15  A  cm~2  and  at  0.3  V  are  shown  in  Fig.  7.  For  both 
the  constant  current  and  constant  voltage  start-up  processes,  the 
ice  distributions  are  similar,  with  slightly  more  ice  on  the  sides,  and 
the  ice  amounts  in  the  different  cells  are  generally  similar  during 
the  whole  start-up  processes,  which  is  also  consistent  with  the 
results  shown  in  Fig.  5.  The  ice  volume  fraction  discrepancies 
among  the  side  cells  and  middle  cells  become  more  notable  during 
the  cold  start  processes,  and  this  agrees  with  the  trend  of  the  stack 
temperature  distribution  shown  in  Fig.  4. 

3.3.  External  heating  and  variable  heating  and  load  control  (VHLC) 

In  order  to  improve  the  cold  start  performance  of  PEMFC  stack, 
external  heating  combined  with  thermal  insulation  is  performed 
and  optimized  based  on  a  20-cell  stack.  The  optimization  is  focused 
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Fig.  6.  Evolutions  of  average  cell  voltage/current  density  and  average  ice  volume 
fraction  in  cathode  CLs  of  different  stacks  starting  from  -20  °C  at  (a)  0.15  A  cm“2  and 
(b)  0.3  V. 
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Fig.  7.  Ice  volume  fractions  in  cathode  CLs  of  different  cells  in  the  50-cell  stack  during 
the  failed  unassisted  cold  start  processes  from  -20  °C  at  (a)  0.15  A  cnr2  and  (b)  0.3  V. 


on  the  control  of  the  heating  location  and  the  load  on  the  active 
cells,  so  called  the  variable  heating  and  load  control  (VHLC).  To 
compare  the  different  cases  consistently,  the  total  external  heating 
power  is  set  as  a  constant  of  100  W  (around  20%  of  the  power  that 
the  20-cell  stack  can  generate  in  normal  operation),  which  can  be 
supplied  by  an  auxiliary  power  source  such  as  a  super  capacitor  or  a 
battery. 

Fig.  8  shows  the  evolutions  of  stack  voltage  and  average  ice 
volume  fraction  in  cathode  CLs.  The  external  heating  is  performed 
uniformly  on  each  individual  cell,  and  for  comparison,  the  cold  start 
processes  with/without  external  heating  and  with/without  thermal 
insulation  on  the  endplate  surfaces  are  all  considered.  With 
external  heating,  the  ice  formation  process  is  slower,  the  stack 
voltage  is  higher,  and  the  cold  start  process  lasts  longer.  The  ther¬ 
mal  insulation  improves  the  cold  start  performance  less  signifi¬ 
cantly  than  the  external  heating.  However,  all  the  cold  start 
processes  shown  in  Fig.  8  are  finally  failed.  Therefore,  for  successful 
cold  start,  either  the  external  heating  power  needs  to  be  increased, 
or  the  heating  method  needs  to  be  optimized. 

A  novel  variable  heating  and  load  control  (VHLC)  method  is 
proposed  in  this  study.  The  main  idea  of  this  method  is  to  only 
apply  load  to  the  still-active  cells  (not  fully  covered  by  ice),  and  to 
simultaneously  apply  external  heating  to  certain  cells  inside  stack. 
This  method  can  be  easily  implemented  by  placing  the  external 
heating  pads  at  different  locations  of  the  stack  and  modifying  the 


external  circuits.  Based  on  this  method,  the  external  heating  power 
is  expected  to  be  used  more  efficiently,  and  the  stack  self-heating 
ability  is  expected  to  be  enhanced.  The  still-active  cells  assisted 
by  the  external  heating  power  may  start-up  successfully  and 
eventually  melt  the  ice  in  the  dead  cells  (cells  fully  blocked  by  ice). 
To  examine  this  method,  the  same  heating  power  of  100  W  is  used 
to  compare  with  the  results  in  Fig.  8. 

By  comparing  with  Fig.  8,  the  cold  start  process  in  Fig.  9  has 
thermal  insulation  on  the  endplate  surfaces,  and  the  external 
heating  is  applied  on  the  two  end  cells  on  both  sides  of  the  stack. 
The  evolutions  of  the  stack  voltage  and  the  ice  volume  fractions  of 
the  three  cells  at  different  locations  of  the  stack  are  shown  in 
Fig.  9(a).  It  can  be  noticed  that  the  stack  successfully  starts  up 
from  -20  °C,  though  a  period  with  relatively  low  stack  voltage.  In 
the  first  50  s,  the  stack  voltage  keeps  increasing,  because  the 
membranes  are  humidified  by  the  product  water  and  the  CLs  are 
not  fully  covered  by  ice  yet.  At  about  50  s,  the  stack  voltage  sharply 
decreases  by  about  7.2  V,  because  most  of  the  cells  failed  at  almost 
the  same  time,  which  can  be  indicated  from  the  ice  volume  fraction 
curves  for  Cell  5  (5  cells  away  from  the  endplate)  and  Cell  10  (in  the 
middle  of  the  stack).  The  stack  voltage  starts  increasing  again  at 
about  60  s,  and  the  step  increment  is  about  0.5  V  for  about  every 
14  s,  indicating  that  the  dead  cells  are  active  again  one  by  one.  For 
the  ice  distribution  in  Fig.  9(b),  it  can  be  observed  that  the  ice 
volume  fractions  in  the  16  middle  cells  are  almost  the  same  and 
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Fig.  8.  Evolutions  of  stack  voltage  and  average  ice  volume  fraction  in  cathode  CLs  of 
the  20-cell  stack  starting  from  -20  °C  at  0.15  A  cm-2  with  the  two  end  surfaces  with/ 
without  thermal  insulation  and  with/without  external  heating  (external  heating  is 
evenly  distributed  on  each  individual  cell). 


become  1  all  at  about  50  s.  On  the  other  hand,  the  4  cells  at  the  two 
ends  of  the  stack  keep  working  due  to  the  external  heating.  For  the 
temperature  distribution  in  Fig.  9(c),  the  temperatures  are  the 
highest  at  the  two  ends  of  the  stack  due  to  the  external  heating. 
However,  it  takes  about  245  s  to  melt  the  ice  in  the  middle  cell. 

By  comparing  with  Fig.  9,  the  change  in  Fig.  10  is  that  the 
external  heating  is  applied  to  the  individual  cell  in  the  middle  of  the 
stack  (Cell  10).  Accordingly,  the  ice  starts  melting  in  Cell  10  first.  It  is 
suggested  that  the  cold  start  process  in  Fig.  10  is  more  optimal  than 
Fig.  9,  because  the  lowest  stack  voltage  during  the  cold  start  process 
in  Fig.  9  is  about  2.2  V,  and  it  is  about  2.8  V  in  Fig.  10,  indicating  that 
there  is  one  more  cell  working. 

It  can  be  noticed  that  the  6  cells  in  the  middle  are  always  active 
during  the  cold  stat  process  in  Fig.  10.  Therefore,  for  the  cold  start 
process  in  Fig.  11,  the  external  heating  is  evenly  distributed  to  the  6 
middle  cells.  The  lowest  stack  voltage  is  about  3.4  V  during  the  cold 
start  process  in  Fig.  11,  which  is  about  0.6  V  higher  than  the  one  in 
Fig.  10,  indicating  that  the  number  of  cells  working  is  increased.  In 
addition,  the  total  time  for  every  cell  to  reach  0  °C  is  reduced  for 
more  than  10  s,  and  the  temperature  is  more  evenly  distributed  for 
the  cold  start  process  in  Fig.  11. 

The  results  in  Figs.  8—11  demonstrate  that  with  the  same 
external  heating  power,  the  novel  variable  heating  and  load  control 
(VHLC)  method  can  improve  the  cold  stat  performance  signifi¬ 
cantly,  which  is  critically  important  for  PEMFC  in  automotive 
applications. 

4.  Conclusion 

In  this  study,  a  cold  start  model  for  proton  exchange  membrane 
fuel  cell  (PEMFC)  stacks  is  developed,  and  a  novel  start-up  method, 
variable  heating  and  load  control  (VHLC),  is  proposed  and  evalu¬ 
ated.  The  main  idea  is  to  only  apply  load  to  the  neighboring  still- 
active  cells,  and  to  apply  external  heating  to  certain  cells  inside 
the  stack  simultaneously  (load  is  not  applied  to  the  cells  fully 
blocked  by  ice,  although  these  cells  can  gain  heat  from  neighboring 
cells).  Simulations  of  the  cold  start  processes  from  -20  °C  are 
performed,  with  the  initial  water  content  in  ionomer  of  2.5  (rep¬ 
resenting  the  well  purged  condition  before  the  cold  start  pro¬ 
cesses),  and  the  initial  ice  volume  fraction  of  0.  The  results  show 
that  the  stacks  with  more  cells  can  reach  higher  temperatures,  and 
for  the  stacks  with  more  than  20  cells,  this  effect  becomes 
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(b)  ice  volume  fractions  in  cathode  CLs  of  different  cells,  and  (c)  temperature  distri- 

surfaces  thermally  insulated  and  with  an  external  heating  power  of  100  W  (external 
heating  is  applied  to  the  two  end  cells  on  both  sides). 
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Fig.  10.  Evolutions  of  (a)  stack  voltage  and  average  ice  volume  fraction  in  cathode  CLs, 
(b)  ice  volume  fractions  in  cathode  CLs  of  different  cells,  and  (c)  temperature  distri¬ 
butions  of  the  20-cell  stack  starting  from  -20  °C  at  0.15  A  cm*2  with  the  two  end 
surfaces  thermally  insulated  and  with  an  external  heating  power  of  100  W  (external 
heating  is  applied  to  the  individual  cell  in  the  middle  of  the  stack). 


Fig.  11.  Evolutions  of  (a)  stack  voltage  and  average  ice  volume  fraction  in  cathode  CLs, 
(b)  ice  volume  fractions  in  cathode  CLs  of  different  cells,  and  (c)  temperature  distri¬ 
butions  of  the  20-cell  stack  starting  from  -20  °C  at  0.15  A  cm*2  with  the  two  end 
surfaces  thermally  insulated  and  with  an  external  heating  power  of  100  W  (external 
heating  is  applied  to  the  6  middle  cells  of  the  stack). 
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insignificant.  With  the  VHLC  method,  the  stack  voltage  first  in¬ 
creases  due  to  the  membrane  hydration,  then  decreases  because 
some  the  cells  are  fully  blocked  by  ice,  and  finally  the  dead  cells  are 
active  again  one  by  one  for  successful  start-up.  With  proper  control 
of  the  VHLC  method,  it  is  demonstrated  that  the  cold  stat  perfor¬ 
mance  can  be  improved  significantly,  which  is  critically  important 
for  PEMFC  in  automotive  applications. 
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